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ABSTRACT 

We present a combined analysis of XMM-Newton, Chandra and Rosat observa- 
tions of the isolated neutron star RXJ0720. 4-3125, spanning a total period of ~ 7 
years. We develop a maximum likelihood periodogramme for our analysis based on 
the AC-statistic and the maximum likelihood method, which are appropriate for the 
treatment of sparse event lists. Our results have been checked a posteriori by folding a 
further BeppoSAX-dataset with the period predicted at the time of that observation: 
the phase is found to be consistent. 

The study of the spin history and the measure of the spin-down rate is of ex- 
treme importance in discriminating between the possible mechanisms suggested for 
the nature of the X-ray emission. The value of P, here measured for the first time, 
is k 10~ 14 s/s. This value can not be explained in terms of torque from a fossil disk. 
When interpreted in terms of dipolar losses, it gives a magnetic field of B ss 10 G, 
making also implausible that the source is accreting from the underdense surround- 
ings. On the other hand, we also find unlikely that the field decayed from a much 
larger value (B « 10 15 G, as expected for a magnetar powered by dissipation of a 
superstrong field) since this scenario predicts a source age of « 10 yrs, too young to 
match the observed X-ray luminosity. The observed properties are more compatible 
with a scenario in which the source is w 10 6 yrs old, and its magnetic field has not 
changed substantially over the lifetime. 

Key words: Stars: neutron; stars: oscillations; pulsars: general; magnetic fields. 



1 INTRODUCTION 

RXJ0720. 4-3125 is a nearby, isolated neutron star (NS) that 
was originally discovered with ROSAT during a systematic 
survey of the galactic plane (Haberl et al., 1997). The source 
exhibits all the common characteristics of the other six 
ROSAT NS candidates (hereafter dim NSs, see e.g. Treves 
et al., 2000 for a review): a blackbody-like, soft spectrum 
with kT ~ 80 eV; an exceedingly large X-ray to optical flux 
ratio; a low X-ray luminosity, Lx ~ 10 30 — 10 31 ergs -1 ; a 
low column density and no evidence for a binary compan- 



* Based on observations obtained with XMM-Newton, an ESA 
science mission with instruments and contributions directly 
funded by ESA Member States and the USA (NASA). 



ion. An optical counterpart has been identified (Motch & 
Haberl, 1998; Kulkarni & van Kerkwijk, 1998). 

Pulsation associated with the spin period has been ob- 
served in four dim NSs: RXJ0720.4-3125 (P » 8.39 s), 
RXJ0420.0-5022 (P « 22.69 s), RXJ0806.4-4123 (P » 
11.37 s) and RBS1223 (P « 5.16 s). The pulse shape of 
RXJ0720.4-3125 is almost sinusoidal (Haberl et al, 1997; 
Cropper et al, 2001). 

The true nature of the mechanism powering the X-ray 
emission from dim NSs is still unclear. Until a few years 
ago, these sources were thought to constitute a stand alone 
class and two mechanisms were proposed for their emission: 
either accretion from the interstellar medium onto an old 
neutron star or release of thermal radiation from a younger, 
cooling object (see e.g. Pavlov et al., 1996; Treves et al., 
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2000; Motch, 2001 and references therein). More recently, it 
has been noted that the period of RXJ0720. 4-3125 is some- 
what unusual and falls in the very narrow range in which the 
anomalous X-ray pulsars (AXPs) periods cluster (P ~ 6—12 
s, see e.g. Mereghetti, 2000). This hints at a possible evo- 
lutionary link between dim NSs, AXPs, and soft gamma- 
ray repeaters (SGRs, see e.g. Heyl & Kulkarni, 1998; Heyl 
& Hernquist, 1998; Colpi et al., 2000; Alpar, 1999, 2000, 
2001). Two kind of "unified" scenarios have been then pro- 
posed. In the first one, the common mechanism powering 
the three classes of objects is dissipation of a decaying, su- 
perstrong magnetic field (B> 10 14 — 10 15 G). In this case 
dim NSs represent the descendants of SGRs and AXPs, and 
RXJ0720. 4-3125 may be one of the closest old magnetars. 
Alternatively, all three classes may contain neutron stars 
with lower (canonical) magnetic field (B w 10 12 G) endowed 
by a fossil disk (Alpar, Ankay & Yazgan, 2001; Alpar, 2001). 
In this interpretation, dim NSs would be in the propellor 
phase and would be the progenitors of AXPs and SGRs, the 
latter having entered an accretion phase. 

In any model involving magnetars, proton cyclotron fea- 
tures are expected to lie in the X-ray range, while electron 
cyclotron lines are expected in the same band for canonical 
magnetic field. According to Zane et al. (2001), for surface 
magnetic fields strengths of ~ 10 14 — 10 15 G spectra exhibit 
a distinctive absorption feature at the proton cyclotron en- 
ergy ~ 0.63(B/10 14 G) keV. The required resolving power 
is w 100, therefore the detection of this feature is well 
within the capabilities of XMM-Newton grating spectrome- 
ters. Recently, Paerels et al. (2001) have presented spectra of 
RXJ0720.4-3125 using XMM-Newton, and found that there 
is no evidence for absorption lines and edges in the X-ray 
spectrum. Unless different atmospheric effects may concur in 
suppressing the feature, when taken straightforward the ab- 
sence of electron or proton cyclotron resonances in the RGS 
band excludes a range of average magnetic field strengths, 
B « (0.3 - 2) x 10 11 G and (0.5 - 2) x 10 14 G. 

Based on the same XMM-Newton observation, Crop- 
per et al. (2001) presented the pulse-shape analysis of 
RXJ0720. 4-3125, modelling spin pulse intensity and hard- 
ness ratio profiles. By assuming that the source of the flux 
variation is the changing visibility of the heated magnetic 
polar caps, they derived an upper limit on the cap size and 
showed that a polar cap larger than ~ 60° — 65° can be 
rejected at a confidence level of 90%. Whatever the mech- 
anism, the X-ray emitting region is therefore confined to 
a relatively small fraction of the star surface. It is worth 
noticing that such a small hot region can not be explained 
only in terms of non-uniformity in the distribution of the 
surface cooling temperature induced by an high magnetic 
field (Greenstein & Hartke 1983; Page 1995). The tempera- 
ture gradient which is induced by the B-dependence of the 
thermal conducivity of the neutron star crust is relatively 
smooth and, even in the (unrealistic) limit case B — > oo 
the associated blackbody luminosity drops only by a factor 
1/2 at ~ 60° and by an order of magnitude at ~ 77°. By 
performing pulse phase spectroscopy, Cropper et al. (2001) 
also found that the observed hardness ratio is softest slightly 
before flux maximum. The same characteristic was later dis- 
covered by Perna et al. (2001) in the spectra of AXPs. Crop- 
per et al. (2001) suggested two possible explanations for this 
effect: their best-fitting model is based on radiation beam- 



ing, while an alternative one assumes a spatially variable 
absorbing matter, co-rotating in the magnetosphere. The 
latter may be indeed the case if the star is propelling matter 
outward (Alpar, Ankay & Yazgan 2001). 

Further information about the nature of this puzzling 
source can be obtained by the spin history. In magnetars 
magneto-dipolar radiation will cause rapid spin-down at a 
rate P « 10~ n (B/10 14 G) 2 /P ss~\ and it has been the pos- 
itive detection of a secular spin-down of this order that has 
suggested the association of AXPs and SRGs with ultra- 
magnetized NSs (Kouveliotou et al., 1998, 1999; Thomp- 
son et al., 2001). The preliminary measure of P published 
by Haberl et al. (1997) for RXJ0720.4-3125 is uncertain 
to a considerably large value, and does not allow spin-up 
and spin-down to be discriminated. Very recently, Ham- 
baryan et al. (2002) presented the first evidence of large 
spin-down rate in a dim NS, RBS1223. Their measure is 
partially based on ROSAT data where the detection of 
the period was not highly significant, and also the value 
(P = 1.35+^7 x 10" 11 s/s) is still compatible with being 
due to torque from a fossil disk. However, it is certainly 
worth noticing that, when taken face value and interpreted 
in terms of dipolar rotational losses, this implies the presence 
of a superstrong magnetic field, B ~ (3.5 — 6.5) x 10 14 G. 

A similar measure of P for the closest pulsating candi- 
date, RXJ0720. 4-3125, is therefore of extreme importance, 
as well as the accurate tracking of its spin history. 

This is the first opportunity to investigate the pulse 
timing in the XMM-Newton data, since at the time of the 
previous papers (Paerels et al. 2001, Cropper et al., 2001) 
the timing correlation files were not sufficiently well de- 
termined. Here we present a combined analysis of XMM- 
Newton, Chandra and Rosat data, spanning a period of ~ 7 
years. 



2 OBSERVATIONS 

The log of observations used for this study is given in Ta- 
ble §. 

RXJ0720.4-3125 has been observed twice by XMM- 
Newton: first on 2000 May 13, during the Calibra- 
tion/Performance Verification phase (Paerels et al. 2001, 
Cropper et al. 2001) and later on 21 November 2001, again 
for calibrations. For the analysis presented in this paper we 
have used data from all three EPIC cameras on both epochs. 
The data were reduced using the XMM-Newton Scientific 
Analysis System (SAS) version 5.1. Source counts in the 
range 0.12 to 1.2 keV (PN) and 0.1 to 2.0 keV (MOS) were 
extracted in an aperture of 40 arcsec for the PN XOOa, 30 
arcsec for PN XOOb, and 30 arcsec for the MOS cameras. 
Events with patterns 0-4 (single and double events) were 
selected for MOS, while patterns 0-12 (all the valid PN pat- 
terns) were retained for PN. Data arriving during episodes 
of high background which were experienced in the 2000 May 
13 observation were excluded. Also in this observation, some 
PN event times were noted not to be integer multiple of the 
frame times, but delayed by 1 second: such events (in to- 
tal about 34 ks from the 50 ks) were corrected in the event 
list. Finally, times were converted to Barycentric Dynamical 
Time (TDB) at solar system barycentre. 

The Chandra observations were made over three days 
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Date 


Observatory 


Instrument 


Exposure 


Exposure 


Effective 


Label 








identification 


duration 


exposure 












(s) 


(s) 




1 qqs Sen 27 


Rosat 


PSPC 


rr>300338n00 


11980 


3221 


R93 


1996 Apr 25 


Rosat 


HRI 


rh300508n00 


7838 


3566 


R96a 


1996 May 7 


Rosat 


HRI 


rhl80100n00 


7743 


3125 


R96b 


1996 Spd 27 


Rosat 


HRI 


rh300508a01 


1498 


1409 


R96c 


1996 Nov 3 


R 


HRI 


rh400884n00 

1 11t:UVJOUt:11UU 


65698 


33569 


R96d 


1998 Sep 27 


Rosat 


HRI 


rh400944n00 


460195 


3566 


R98 


2000 Feb 1 


Chandra 


HRC-S(LETG 1st order) 


348+349+745 


305528 


37635 


ChOO 


2000 May 13 


XMM-Newton 


MOS1 + thin filter 


0124100101-001 


61648 


61648 


XOOa 






MOS2 + thin filter 


0124100101-002 


61648 


61648 








PN + thin filter 


0124100101-003 


62425 


62425 




2000 Nov 21 


XMM-Newton 


MOS1 + medium filter 


0132520301-007 


17997 


17997 


XOOb 






MOS2 + medium filter 


0132520301-008 


17994 


17994 








PN + medium filter 


0132520301-003 


25651 


25651 




1997 Mar 16 


BcppoSAX 


LECS 


LECS.20079001 


99418 




S97 



Table X. The ROSAT, Chandra and XMM-Newton observations of RXJ0720.4-3125 used in this paper. 



starting on 2000 Feb 1 using the HRC-S with LETG. The 
events were extracted from zero-order images of the source 
(namely, from circles of 2 arcsec radius centered at the star's 
position). The dispersed data (I s * and higher orders) are 
strongly contaminated by background (about 50-60%) and 
are not useful for this task. On the other hand, the back- 
ground contamination in the zero-order data is negligible, 
about 2-3%. Also, the zero-order data are not piled-up: this 
effect can appear only at larger count rates (~ 10 5 ct/s). 

The Rosat observations were taken over a period of five 
years in 1993, 1996 and 1998. An earlier report of the analy- 
sis of the 1993 and 1996 data is in Haberl et al. (1997). The 
data were re-extracted and re-analysed using the EXSAS 
software system to make use of the most recent knowledge 
of spacecraft clock corrections. The 1998 data were taken 
shortly before a spacecraft emergency safe-hold, and do not 
benefit from a clock calibration: timing corrections are there- 
fore extrapolated from the last available calibration. Tim- 
ings may therefore be incorrect up to several seconds (W. 
Becker, private communication): this means that the refer- 
ence of this observation to the rest of the dataset is prob- 
lematic. However, no significant drifts are expected over the 
duration of the observation itself, so that the data are still 
useful in a standalone analysis. 

Although the countrate is too low for our main analysis, 
we have also extracted the SAX LECS observation of RX 
J0720. 4-3125 from the SAX data archive (the source was 
not visible in the other SAX instruments) for use as an a 
posteriori check on our derived timing parameters. Here we 
used the event list from the pipeline, correcting the times to 
the solar system barycentre using the FTOOL earth2sun. 
This neglects the light travel time from Earth centre to the 
satellite in low Earth orbit, ~ 20 ms, which is negligible for 
the purposes of our a postiori check. 

The major datasets in Table ^ are from the two XMM 
Newton observations which have high count-rates (particu- 
larly the EPIC-PN), and the long 1996 Nov 3 Rosat data. 
The 1998 Rosat and the Chandra observations are valuable 
by nature of their several day durations. 
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Figure 1. The Rayleigh transform of the XOOa PN dataset over 
the period range 100 s to 2 s. This shows only the single frequency 
at 8.391 s (0.11917 Hz), as reported in Haberl et al. (1997). 

3 TIME SERIES ANALYSIS 

Our data originate from instrumentation with widely differ- 
ing sensitivities: typical count rates vary from 1 count every 
~ 3 s for Rosat HRI to ~ 6 counts/s for XMM-Newton PN. 
However, none of these count rates is sufficiently high for 
a normal distribution of counts to be expected in a time 
bin of adequate time resolution. The body of literature ad- 
dressing discrete Fourier Transforms (for example Deeming 
1975, Scargle 1982, Schwarzenberg-Czerny 1998) is not di- 
rectly applicable for our analysis. The treatment of sparse 
data and event list data is generally accomplished using 
the Rayleigh Transform (for example de Jager 1991, Mardia 
1972). It should be noted that this transform is numerically 
the same as the discrete Fourier Transform. 

It will be important in our analysis to define confidence 
intervals to the derived quantities, in particular the period. 
In the case where sufficient counts per sample are available, 
this is generally achieved by least squares fits to the data 
(for example Kurtz 1985), but this is not possible for sparse 
or event list data. In this case it is necessary to use the more 
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Figure 2. (Left) Maximum likelihood periodogrammes (MLP) for three long datasets, R96d, XOOa (PN) and XOOb (PN), showing the 
unambiguous detection of a periodicity at 8.391 s. These constrain the selection of the strongest and second-strongest dips in the MLPs 
for the R98 and ChOO datasets respectively (right). The vertical line denotes a period of 8.39113 s. The 68% and 90% confidence levels 
are at x 2 = 1.0 and 2.71 for one degree of freedom. 



general maximum likelihood technique, which makes no as- 
sumptions on the data distribution (Cash 1979, Lampton, 
Margon & Bowyer 1976). Using this, Cash (1979) defined the 
C- and AC-statistics. The AC-statistic (Cash 1979) has the 
considerable advantage of being distributed as x 2 , s ° that 
the required confidence intervals can be identified directly 
from the plot. 

The maximum likelihood technique depends on a model 
probability distribution, which must be specified a priori. 
Haberl et al. (1997) and Cropper et al. (2001) find the flux 
variation of RXJ0720. 4-3125 to be almost sinusoidal, so in 
our particular case we can take a sinusoidal variation su- 
perimposed on a mean value as an appropriate model. The 
single frequency and absence of harmonics in the Rayleigh 
transform of the XOOa PN dataset shown in Figure hi pro- 
vides the justification explicitly for this approach. The max- 
imum likelihood is then computed simply as equation (35) of 
Bai (1992). The normalisation of the model is unimportant, 
as it results in a constant term in the C-statistic, which then 
is eliminated in the AC-statistic. In our case we have two 
free parameters: normalised amplitude and phase. In gen- 
eral, when performing maximum likelihood fits, the best fit 
free parameters have to be determined by whatever method 
may be most appropriate. However, in the case of sinusoids, 
as noted by Bai (1992), their orthogonality allows them to be 
determined directly from the Rayleigh Transform (or Lomb- 
Scargle periodogramme). 

The computation of the AC-statistic is therefore only 



slightly more complicated than the periodogramme. It starts 
in the same way, by a scan through period, calculating the 
Rayleigh Transform to provide the amplitude and phase at 
each period. Then at each period the maximum likelihood 
is calculated using equation (35) of Bai (1992). Finally from 
equations (4) and (7) of Cash (1979), AC is simply twice the 
difference between the maximum likelihood at a particular 
period, and the maximum value over the whole frequency 
range. This results in an inverted "periodogramme", with 
minimum value of zero. The permitted period range is that 
between the confidence levels appropriate for \ 2 with the 
appropriate number of degrees of freedom (1 or 2 depending 
on whether P is included as a parameter - see section ^J) . 

Details of the procedure are given in Appendix A. 

With a period of 8.391 s (Haberl et al. 1997), more than 
10 7 cycles have elapsed over the time span of our data set. It 
is therefore essential to incorporate a P term in our analy- 
sis techniques when combining datasets at different epochs. 
The search then becomes a (Po , P) search, where Po is the 
period at the start of the first dataset (the 1993 Rosat ob- 
servations), and the instantaneous period P = Po + P(t — to) 
is computed at each point in the grid. Here to is the time at 
the start of the R93 run. 

As a more general comment, we recommend this max- 
imum likelihood periodogramme (hereafter MLP) for time 
series analysis. The uncertainty in a period determination 
is not evident directly from the width of a peak in a power 
spectrum computed by classical discrete Fourier transform, 



© 0000 RAS, MNRAS 000, 000-000 



Timing analysis of the isolated neutron star RX J0720. 4-3125 5 



Lomb-Scargle periodogramme or Rayleigh transform, as this 
is set by the window function (c.f. Scargie 1982 appendix D). 
However, it is directly available from the MLP AC statistic: 
the inverted peak from a sinusoidal signal in data with high 
signal-to-noise ratio will be narrower than the peak from 
low signal-to-noise data and the \ 2 can be read directly 
from the vertical axis (see e.g. Figure ^). Furthermore, as 
noted above, the MLP technique requires an appropriate 
model, which may be unknown a priori. In this case the 
signal can be reconstructed from the harmonic content of 
the Rayleigh transform or Lomb-Scargle periodogramme, or 
from the phase folded data, and an appropriate model gen- 
erated. This model can then be used in the calculation of 
the MLP, incorporating all of the harmonic information ex- 
plicitly. 
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R98 


2450923.523982 
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0.000013 


ChOO 


2451575.771501 
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0.000020 


XOOa PN 


2451677.604110 


8.391133 


0.000021 


XOOa MOS1 


2451677.614469 


8.391103 


0.000043 


XOOa MOS2 


2451677.614474 


8.391120 


0.000040 


XOOb PN 


2451870.308664 


8.391181 


0.000050 


XOOb MOS1 


2451870.390808 


8.391041 


0.000167 


XOOb MOS2 


2451870.390823 


8.391032 


0.000181 



4 MEASURE OF THE SPIN PERIOD 

The only measure of P and P presented in the literature 
for RXJ0720.4-3125 was given by Haberl et al. (1997) using 
Rosat data. The period was determined to be 8.39115 ± 
0.00002 seconds from the 1996 Nov 3 dataset. Based on a 
linear fit of the data spanning 3 years, these authors derived 
a 90 % confidence range of -6.0 x 10~ 12 to 0.8 x 10~ 12 for 
P, compatible with no period change. 

It is computationally unfeasible to explore this range 
in the PoP plane in a single analysis of the entire dataset 
to determine refined values for Po and P. A reduction in 
the range is required. This can be accomplished initially by 
extending the Haberl et al. (1997) analysis to include the 
AMM-JVewton data. This is sufficiently accurate to identify 
the correct peak in the alias patterns in the extended du- 
ration R98 and ChOO data, leading to a further limiting of 
the permitted range in the PP plane. The next stage is a 
grid search in P and P over this restricted range for the 
1993 and 1996 combined Rosat data only. The final stage is 
a grid search using all data in the restricted zone for which 
there are acceptable solutions from this Rosat analysis. 



4.1 Individual datasets 

We begin by performing an MLP assuming P = on each of 
the longer pointing (R93, R96d, XOOa, XOOb) in the period 
range (8.375-8.400) s, incorporating the value of 8.39115 s 
found by Haberl et al. (1997). The results are shown in Fig- 
ure | (left). As we can see, there is no ambiguity in the period 
determinations from the R96d, XOOa and XOOb datasets; the 
best fit frequencies and their uncertainties are reported in 
Table ^. A linear least square fit using the 68% formal errors 
from the MLP results in a value of P = 8.39113±0.00011 s, 
P — 0.0 ± 5.5 x 10~ 13 s/s, where, as throughout, Po is ref- 
erenced to the start of the R93 run. 

The proximity in time between the ChOO and XOOa ob- 
servations, together with the above upper limit in P, per- 
mits an unambiguous determination of the second high- 
est peak in the ChOO power spectrum (Figure ^, right). 
The same is true for the R98 observations, in which the 
highest peak is selected. Adding these to the linear least 
squares fit results in a value of Po = 8.39107 ± 0.00005 and 
P = 2.7 x 10~ 13 ± 2.5 x 10~ 13 . As may be expected, the P 



Table 2. Best fit periods to the longer individual datasets, as- 
suming P = 0. 



and P are highly correlated. The 68, 90 and 99% confidence 
intervals are shown in Figure |j| 

The duration and high count-rate of the XOOa PN run 
alone provides a strong constraint in the (Po , P) plane. The 
MLP 68% and 90% confidence intervals are also shown in 
Figure ^. This indicates that acceptable values for Po and 
P lie within the range 8.39106 to 8.39115 s and -1 to 
+2 x 10~ 13 s/s, respectively. We checked the effect of a 10 -7 
error in the XMM-Newton clock calibration, as reported re- 
cently (Kuster, M., 2001, reported at "New Visions of the 
X-ray Universe in the XMM-Newton and Chandra Era"). 
The change in the confidence region is negligible. 

4.2 Combined Rosat dataset 

The Po P range identified in Figure |§] is sufficiently restricted 
for an MLP to be performed on the combined 1993 and 
1996 Rosat datasets. The use of combined datasets requires 
coherent phasing to be maintained over the whole dataset 
(coherent phasing was required only within each of the in- 
dividual datasets used until this point). The (Po ,P) plane 
was searched over between 8.39095 < Po < 8.39120 s, and 
-4 x 10~ 13 < P < 8 x 10~ 13 s/s, exceeding the parame- 
ter range for the 90% confidence limits from the linear least 
squares fit. The resulting 68% and 90% confidence contours 
break up into small regions (aliases) distributed along lines 
in the (Po , P) plane. These are also shown in Figure |3j 

As is evident, there is overlap between the Rosat 90% 
confidence interval contours and the 99% contours of the 
linear least squares fit to the individual data (but not quite 
at the 90% level) . This overlap region is also consistent with 
the XOOa PN confidence intervals. 

4.3 Full dataset 

With the further restriction available from the fit to the 
Rosat 1993 and 1996 data subset in Figure |§| a coherent 
MLP search was made to all of the data (excluding the 
uncertain R98 run). The search was made for 8.391050 < 
Po < 8.391150 s and -3 x 10" 13 < P < 3 x 10" 13 s/s 
within the parallel lines Po = aP + b defined by intercepts 
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Figure 3. The 68%, 90% and 99% contours for Po and P linear least squares fit to the individual best fit periods in Table || (continuous 
elliptical regions). Overplotted on the same scale are the 68% and 90% contours from the MLPs to the XOOa PN (defined by continuous 
parallel lines) and the combined 1993 and 1996 Rosat datasets. The latter are broken up into tiny elliptical regions corresponding to the 
68% and 90% confidence levels for the different aliases: see the enlargement inset in the top right hand part of the figure. 



b = 8.391097 s and 8.391105 s and slope a = -1.75 x 10 8 
s, encompassing the Rosat 1993/1996 contours. Care was 
taken to ensure adequate sampling of the (Po , P) plane. 
This identified two pairs of values, (Po , P) which cannot 
be discriminated between on statistical grounds. These are 
(Po , P) = 8.39109273 s, 5.409 x 10~ 14 s/s and 8.39109148 s, 
3.749 x 10" 14 s/s (Table |). All other periods can be ex- 
cluded at the 95.4% level, but two other possibilites cannot 
be excluded at the 99% level: these are also given in Table 
|§] for completeness. The contours for the two best fits are 
given in Figure ^ 

We have folded the data on both PoP solutions (1) and 
(2) in Table ^| to check that the relative phasing of all in- 
dividual runs is correct. We show the folded data on these 
periods in Figure []. The data from runs R96a and R96b have 
been combined to increase the signal-to-noise ratio, and we 
have not shown the data from R98 because of the uncer- 
tainty in the time reference for this run (section 2). 



Label 


Po 


P 




Ax 2 




(s) 


(s/s) 






(1) 


8.39109273 


5.409 x 10" 


14 




(2) 


8.39109148 


3.749 x 10- 


14 


1.3 


(3) 


8.39108624 


6.902 x 10- 


14 


6.2 


(4) 


8.39108748 


8.564 x 10" 


14 


6.3 



Table 3. The four best fit (P , P) values. Refer to Figure | for 
the confidence interval contours. Ax 2 is the difference between 
the x 2 of a given solution, and that of solution (1). Fits (3) and 
(4) are formally excluded at the 95.4% level. 



Now performing the a postiori check with the BeppoSax 
data, we find that these data phase correctly with the main 
datasets for PqP (1), but not (2). This is shown in the lowest 
panel of Figure in each case. This allows us to select PqP 
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(i) 



8.39109280 - 



8.391 3>)2~j - 



8.39109270 - 



8.39109265 - 



8.39109260 - 




(2) 



8.39109160 - 



8.39'09155 - 



i.39109145 - 



8.39109140 - 




5.38 5.40 5.42 

Pdot x10~'* (s/s) 



3.70 3.72 3.74 3.76 
Pdot x10~" (s/s) 



Figure 4. The 68% and 90% MLP contours for Pn and P to the complete dataset except for the R98 data, for the two acceptable 
(Pq i P) values (1) and (2). Note that the contours in (2) are referenced to the best global fit derived for (1). 



(1) as the most likely timing parameters for RX J0720.4- 
3125. In any case, both acceptable fits to the data have 3 x 
lCT 14 < P < 6 x 10" 14 . For the purposes of our further 
discussion, the difference between the P in (1) and (2) is 
not important. 



5 DISCUSSION 

The refined value of P reported in this paper is consistent 
with, but two orders of magnitude lower than the extrema 
of the range reported by Haberl et al. (1997). The first im- 
plication is that RXJ0720. 4-3125 is unlikely to be spinning 
down under propeller torques. If the propeller is due to the 
torque exerted by the interstellar medium, it would be (e.g. 
Colpi et al., 1998) 



-8 9/13 

n ' v 



27/13 R 8/13 p 21/13 
10 -"12 " 



S 



(1) 



where n is the external density in cm~ 3 (n ~ 1 for the in- 
terstellar medium), B12 = B/(10 12 G) and via is the star's 
velocity normalized to 10 km/s. By using the values of P 
and P of RXJ0720.4-3125, we obtain B12 « 16n~ 9/8 -u 2 78 . 
Propeller spin-down dominates over dipolar losses if B12 < 
25v / w« 10 3 ^ 2 and this, combined with the previous equation, 
constrains the star's velocity to extremely small values: 
V10 Si n 1 / 3 . Moreover, this scenario makes the matching of 
thermal and dynamical time scales more difficult. On the 
other hand, the fact that the propeller contribution is neg- 
ligeable supports larger values of the star's velocity. 



In the "unified" model suggested by Alpar (2001), the 
observed spin-down is associated with propeller effect by a 
fossil disk. In this case, under the assumption that the X-ray 
luminosity of the source is supplied by energy dissipation 
(frictional torque), upper and lower bounds on £1 can be 
derived by eq.(2) of Alpar (2001). For RXJ0720.4-3125 the 
observed luminosity is L w 2x 10 31 (d/100 pc) 2 erg cm' 2 s _1 , 
where dioo = d/100 pc and d is the distance (Haberl et al. 
1997). This gives: 



2 x 10~ 12 d 2 nn ^ tl< 2 x 10~ 10 d 2 



rad 



that in turn translates in P — (P 2 /2vr)n 



2 x 10 n d? 00 ~ P~ 2 



x 10" 



J 2 s 
3 100 _ 



(2) 



(3) 



While this scenario is still consistent with the spin-down 
measured for RBS 1223 (Hambarayan et al., 2002), the value 
of P reported here for RXJ0720.4-3125 is well below this 
range. Thus, for RXJ0720. 4-3125, the associated energy dis- 
sipation cannot, alone, account for the source luminosity 
(unless the source is at ~ 5 pc, which is unrealistic). That 
also make less plausible an interpretation of the hardness 
ratio profile in terms of a spatially variable absorbing mat- 
ter, co-rotating in the magnetosphere (Cropper et al. 2001). 
The observed behaviour is more probably explained by the 
angle-dependent properties of the emitted radiation. 

On the other hand, the slow spin-down rate of 
RXJ0720.4-3125 is still considerable. The other plausible 
mechanism which may account for such large and stable val- 
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Figure 5. The datasets folded on P P (1) (left) and (2) (right) in Table |[ 
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ues of P is rotational loss by emission of magnetic radiation. 
For a dipolar magnetic field, 

P » W- 15 B 2 12 /P J , (4) 

which gives a present value of the magnetic field of B = 
2.f3 x 10 13 G Pn. Despite this value being below those typ- 
ically quoted for magnetars (B w 10 14 — 10 15 G), it is still 
extreme and close to the critical value Bq — 4.41 x 10 13 G at 
which quantizing effects start to be important in shaping the 
atmospheric emission. Furthermore, a scenario in which this 
source is powered by accretion from the interstellar medium 
must be ruled out: for the present values of P and B the 
corotating magnetosphere will prevent the incoming mate- 
rial to penetrate below the Alfven radius. 
The corresponding spin-down age is 

t sd = PI (2P) ~ 2.48 x 10 6 yr, (5) 

which, given the numerous uncertainties, is marginally com- 
patible with that inferred by the cooling curves: a few 10 
yrs for a surface temperature of ~ 80 eV (see e.g. Kaminker, 
Haensel & Yakovlev, 2001, Kaminker, Yakovlev & Gnedin, 
2001; Schaab et al. 1997, 1999). The discrepancy is less sig- 
nificant if we notice that: a) the cooling curves are strongly 
dependent on the neutron star mass (Kaminker, Haensel & 
Yakovlev, 2001), and b) what we are probably observing in 
X-rays is a region of limited size which is kept hotter than 
the average star surface, as inferred by the analysis of the 
pulse-shape (Cropper et al., 2001). 

On the other hand, t s d is representative of the true age 
of the source only in the case in which the magnetic field 
remained almost constant during the star evolution. The 
same condition applies for the validity of the cooling curves 
mentioned above, which do not include the extra input of 
energy released in the neutron star in case of field decay. It 
is therefore of fundamental importance to address the field 
evolution: the strength of B and its variation during the 
neutron star history determine, in fact, the actual luminos- 
ity, the lifetime and even the nature of the energy loss from 
the star. A related question is where the field is anchored: 
in the core, penetrating the whole star or confined in the 
crust. There are three meachanisms which are typically pro- 
posed for inducing field-decay: ambipolar diffusion in the 
solenoidal or irrotational mode and Hall cascade (Goldreich 
& Reisenegger 1992, Heyl & Kulkarni 1998, Colpi, Geppert 
& Page 2000) . In reality all the three processes co-exist with 
different timescales, and each of them may dominate the 
evolution depending on B, L at any given time. In absence 
of more detailed computations, Heyl & Kulkarni (1998) and 
Colpi, Geppert & Page (2000) tentatively isolated the three 
processes and proposed some simple, phenomenological rules 
to mimic the evolution in the three regimes. We stress that 
these descriptions are far from being comprehensive of all 
the effects that can modify substantially the results: they are 
used in a first approximation and have the main advantage 

t Note that a twisted magnetosphere may lead to a reduction 
up to an order of magnitude in the inferred polar value of the 
magnetic field. The considerations about the age, however, remain 
unchanged (Thompson, Lyutikov & Kulkarni, 2001). 
' Here and in the following we specify the discussion to solution 
(1) of Tabic | 



B-Decay Mechanism 


Bo 


age 




10 13 G 


(years) 


Hall Cascade 


119.2 


4.5 x 10 4 


Ambipolar diffusion, irrotational mode 


1.9 


3.3 x 10 6 


Ambipolar diffusion, solenoidal mode 


4.1 


1.6 x 10 6 



Table 4. Predicted source age and primordial field for three dif- 
ferent mechanisms of decay, simulated as in Colpi et al. (2000). 
The present values of P and P are those of solution (1) in table 
a. In all cases, the source is assumed to be born with P = 1 ms. 



of having a simple analytical form. By using the expressions 
of Colpi et al. (2000), we have estimated the source age and 
the value of the magnetic field at the birth of the neutron 
star, Bo. Results are shown in Table In all cases the star 
is assumed to be born with a period of 1 ms: results are not 
strongly dependent on this exact value, provided it is far less 
than the present period. 

As we can see, allowing for a mechanism involving very 
fast decay, such as the Hall cascade, we find that the source 
is now ~ 4 x 10 4 yr old, and is born with a superstrong field 
-Bo ~ 10 G. Such a young age is only marginally compat- 
ible with the absence of a remnant and, more important, 
is not compatible with the low observed X-ray luminosity 
(as we compared using not only the standard cooling curves 
mentioned above, but also some cooling curves kindly pro- 
vided by Geppert & Colpi, private communication; these 
latter are computed allowing for the extra-heating due to B- 
decay from Hall cascade and predict larger luminosities than 
the former, making the discrepancy even higher). Underlu- 
minous models have been recently presented by Kaminker, 
Haensel & Yakovlev (2001), who accounted for the enhanced 
neutrino cooling in presence of strong neutron superfluid- 
ity. These solutions may match an age of ~ 10 4 yrs for 
RXJ0720. 4-3125, but, as discussed by the same authors, 
they must probably be rejected since they fail in the com- 
parison with observational data of a sample of other neutron 
stars. 

On the other hand, both mechanisms involving ambipo- 
lar diffusion predict a magnetic field quite stable over the 
source lifetime and close to the actual value. Accordingly, 
the predicted age is ~ 10 6 years in all cases, close to t s d- Be- 
low ~ 10 14 G the cooling curves are not significantly influ- 
enced by decay through ambipolar diffusion (Heil & Kulka- 
rni, 1998), thus, as in the case of constant B discussed above, 
the scenario is compatible with the observed luminosity. The 
larger age is also compatible with the absence of a remnant. 

If our conclusions are valid, the connection between dim 
INS and AXPs is not so obvious. RXJ0720. 4-3125 has a 
strong, but not superstrong, field which is compatible with 
those of the canonical radio-pulsars which have passed the 
death line. On the other hand, having excluded accretion, 
what mechanism causes an X-ray emission concentrated in 
a fraction of ~ 60% of the star surface remains a mystery, 
as well as the related question about the validity of using 
the observed blackbody temperature to locate the source in 
the cooling diagram. The variation of the surface cooling 
temperature with the latitude predicted so far for strong 
fields is not large enough and more complicated explanations 
are required. 
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APPENDIX A: DETAILS OF THE MAXIMUM 
LIKELIHOOD PERIODOGRAMME 

The C-statistic is derived from the maximum likelihood ra- 
tio and is 



C = 2(£- 



E 



rii ln/i) 



(Al) 



(Cash 1979), where E is the sum of the expected number 
of counts according to the model distribution, rii is the ob- 
served number of counts in an interval and Ii is the expected 
number of counts at the time of event ti according to the 
model distribution. N is the total number of counts. For 
event list data, the sampling duration tends to zero and 
rii = 1- The AC-statistic is 



AC — i^Cmiri) p — q (CmaTi)p 



(A2) 



where p is the total number of parameters, and q is the 
number of free parameters. This is distributed as x 2 with 
q degrees of freedom (Cash 1979). (C m in) P is the minimum 
C-statistic when all parameters are varied - the minimum 
over all points in the search grid in the free parameter 
space. (C m in)p-q is the minimum C-statistic at any par- 
ticular point in the search grid. Substituting, 



AC = 2 



rii In h 



22 Ui m Ii 



particular 



.(A3) 



The expected number of counts U in equation (Al) is 
calculated from the model distribution. In our case of a pure 
sinusoidal variation, because the expected arrival of an event 
is directly proportional to the model prediction for time ti, 
we have 



= a (l + Acos(27rfi/P + O ) 



(A4) 



where A is the normalised amplitude, P is the period and 0o 
is the phase. The scaling of this model ao is unimportant, as 
constant factors are eliminated when substituting in equa- 
tion (A3). This model has parameters (A,P, 0o) so p = 3. 
The periodogramme will scan in P, with best fit values for 
A and 0o, so q = 1. In the case where we allow a period 
change, 



P = Po + PU 



(A5) 



where Po is the period at ti = and P is the period change, 
the periodogramme will consist of a scan in Po and P. In 
this case p = 4 and q — 2. 

In the general case of a maximum likelihood statistic, 
an optimisation search is required to obtain the best fit for 
the p — q fitted parameters. Here, in the case of a sinusoid, 
these can be derived directly from the Rayleigh transform, 
as pointed out by Bai (1992), 



1 

N 



V m cos 27rfi / P 




sin 2ivti I P 



(A6) 



where z is the Rayleigh power, and again, rii 
list data. The amplitude is then 



1 for event 
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8.3910 
Period (s) 



8.3915 



8.3920 



Figure 6. The discrete Fourier transform (top) and maximum likelihood periodogramme (bottom) for the XOOa PN and XOOb PN data 
over a restricted period range around the 8.391 s period in RXJ0720. 4-3125. The 68% and 90% confidence levels are indicated at x 2 = 1-0 
and 2.71 for one degree of freedom. 



.4 



and the phase is 



(A7) 



= arctan 



JV 

E 



rii sm2nti/P 




i COS 2nti I P 



This is equivalent numerically to the discrete Fourier trans- 
form (DFT) (Deeming (1975), as can be ascertained by ref- 
erence to Kurtz (1985), equation (A2). 

The procedure is then for each P (which may be cal- 
culated from equation (A5) in a Po ,P search) to calcu- 
late the amplitude and phase through equations (A7) and 
(A8) (which amounts to calculating the Rayleigh or discrete 
Fourier transform), then calculate AC through equations 
(A4) and (A3). As the P grid is searched, the maximum 
m In Ii is stored: AC is then simply twice the differ- 
ence between the calculated value at each period, and this 
recorded maximum. This difference is distributed as \ 2 f° r 
one or two degrees of freedom, depending on whether P 
is used directly, or calculated in a (Po , P) search through 
equation (A5). 

An illustration of the power of the MLP is given in 
Figure |. Here we have taken the XOOa PN and XOOb PN 
data, and computed a discrete Fourier transform in the nar- 
row period range 8.390 to 8.392 s (a much expanded scale 
by comparison with Figure ^), assuming no P. This (upper 
plot) shows a multiplicity of aliases created by the long gap 
in the data, superimposed on the broad peak corresponding 
to the duration of the longest observation (XOOa). These are 



barely resolved in the plot. The lower plot is the MLP with 
the 68% and 90% confidence levels indicated. This eliminates 
most of the aliases, as the narrowness of the distribution is 
related to the precision with which the best fit period can 
(Agpe identified within the broad peak in the upper plot. In the 
case of the high signal-to-noise ratio XOOa and XOOb data, 
this is significantly narrower than the width of the peak, 
which is set by the window function (Scargle 1982) . 
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